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We have constructed maximally-localized Wannier functions for prototype structures of solid molec- 
ular hydrogen under pressure, starting from LDA and tight-binding Bloch wave functions. Each oc- 
cupied Wannier function can be associated with two paired protons, defining a "Wannier molecule" . 
The sum of the dipole moments of these "molecules" always gives the correct macroscopic polariza- 
tion, even under strong compression, when the overlap between nearby Wannier functions becomes 
significant. We find that at megabar pressures the contributions to the dipoles arising from the 
overlapping tails of the Wannier functions is very large. The strong vibron infrared absorption 
experimentally observed in phase III, above ~150 GPa, is analyzed in terms of the vibron-induced 
fluctuations of the Wannier dipoles. We decompose these fluctuations into "static" and "dynami- 
cal" contributions, and find that at such high densities the latter term, which increases much more 
steeply with pressure, is dominant. 
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I. INTRODUCTION 



A. Wannier functions 



The electronic structure of periodic solids is usually 
described, in the independent-electron approximation, in 
terms of the extended Bloch eigenfunctions. An alterna- 
tive representation is provided by the Wannier functions 
(WFs) which are localized, with a typical spread 

of the order of the atomic dimensions; they can be ob- 
tained via a unitary transformation of the Bloch states 
belonging to an isolated band jl] or to a composite group 
of bands |^,|| (i.e., bands that may be connected among 
themselves by degeneracies, but are separated from all 
others by energy gaps). For some purposes the latter 
description is advantageous: for instance, the WFs con- 
structed from the states in the valence bands provide 
an intuitive, "chemical-like" localized picture of bonding 
and dielectric properties of insulators |3|. 

The major drawback of the Wannier representation is 
the strong nonuniqueness of the WFs: their average lo- 
cation, shape, and spread all depend on the arbitrary 
choice of gauge [^,^. In practice, this indeterminacy can 
be resolved by working with the set of WFs which is most 
localized according to some sensible criterion. A certain 
degree of arbitrariness still remains regarding which mea- 
sure of localization to use, and in fact several alternatives 
have been proposed in the literature. We follow the ap- 
proach of Ref. 1^ , which amounts to minimizing the sum 
of the quadratic spreads of the WFs (see Sec. ||). We will 
use density-functional theory in the local density approx- 
imation (LDA) , complemented by a tight-binding analy- 
sis, to investigate in terms of well-localized WFs the elec- 



tronic structure and dielectric properties of compressed 
molecular hydrogen H|. 



B. Compressed molecular hydrogen 

Solid hydrogen under pressure has attracted consider- 
able attention over the years [H, and the main goal has 
been to try to metallize it; this is expected to occur at 
high enough pressures, either by band gap closure in a 
molecular phase, or by molecular dissociation, whichever 
occurs first However, up to the highest pressures 

reached so far 340 GPa), hydrogen appears to remain 
both molecular and insulating Nevertheless, a rich 
phase diagram has emerged, with three distinct phases 
unambiguously identified using Raman and infrared (IR) 
spectroscopy ||]. 

The precise crystal structure of the high-pressure 
phases (phases II and III) has not been determined ex- 
perimentally, and conclusive theoretical predictions have 
proven quite difficult, due to the quantum effects associ- 
ated with the protons. The purpose of the present work 
is not to propose new candidate structures, but rather to 
make some very general points, illustrated on a couple of 
particularly simple prototype structures (Fig. which 
were chosen mainly for clarity. It is hoped that, even 
if none of them turns out to be the correct structure of 
phase III (which is likely to be the case), they manage 
to capture some of its relevant features. For instance, 
ab-initio calculations at megabar pressures and low tem- 
peratures tend to favor structures in which the centers 
of the molecules form an hep or, more generally, a trian- 
gular lattice (possibly with a small distortion [pl]|) 
|l2| . It is believed that in phase III the molecules are ori- 
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FIG. 1. The Cmc2i structure viewed along the c-axis 
(left) and in the yz plane (right). The centers of the molecules 
lie on hep sites, and the molecules in the two sublattices are 
tilted away from the c-axis by opposite angles 9 and —9. The 
C2/m structure is identical except that the two molecules in 
the primitive cell are tilted in the same direction by an angle 
9. The arrows on the right side indicate the directions of the 
dipoles of the two "Wannier molecules" for = 1.52. 



sures, in order to convert from to pressure we use an 
experimental equation of state extrapolated to high pres- 
sures The LDA calculations were performed using 
a plane-wave cutoff of 90 Ry and the bare Coulomb po- 
tential of the protons. The self-consistent calculations 
with 4 atoms per cell used a (11,11,11) Monkhorst- 
Pack mesh for the Brillouin zone sampling. After self- 
consistency was achieved, the "maxloc" WFs were de- 
termined starting from the Bloch states again calculated 
on a (11,11,11) mesh. In all the calculations we have 
used the following parameters for the Cmc2i structure 
described in Fig. |^: rbond — 1.445 a.u., c/a = 1.576, 
and the tilt angle 6 — 54.0°. For C2/m the parameters 
are rbond = 1.456 a.u., c/a — 1.588, and 6 = 69.5°. In 
both cases the structures were obtained by minimizing 
the enthalpy at a fixed LDA pressure of 100 GPa, with 
a resulting density of = 1.52 (which experimentally 
corresponds to about 115 GPa, according to Ref. |p8[). 
The same parameters were used at all other densities. 



entationally ordered, with their axes tilted away from the 
c-axis, as such canted structures tend to be more stable 
|f7|,p| jl0|jr^ . Moreover, the resulting lowering of symmetry 
gives rise to IR-active vibron modes ]l^ , p^ ; indeed, one 
of the signatures of phase III, above 150 GPa, is a strong 
IR absorption peak in the vibron frequency range Jl6| , pTt , 
which contrasts with the much weaker absorption found 
in phases I and II. 



II. MAXIMALLY-LOCALIZED WANNIER 
FUNCTIONS 

A set of WFs {wn(r — R)}, each labeled by a differ- 
ent Bravais lattice vector R, can be constructed from 
the Bloch eigenstates {ipnk} in band n using the unitary 
transformation 



C. Organization 

The paper is organized as follows. In Sec. || we briefly 
review the method used for constructing well-localized 
WFs; in Sec. [II we investigate the permanent dipole 
moments of the "Wannier molecules" as a function of 
crystal structure and pressure; the results are presented 
in terms of a "static effective charge" vector associated 
with each molecule. In Sec. IV we look at the vibron- 
induced fluctuations of those dipoles, which can be quan- 
tified in terms of a "vibron-induced effective charge" vec- 
tor on each molecule. From the effective charges the 
strength of the vibron IR absorption is calculated and 
compared with experiment, and the relative importance 
of the "static" and "dynamical" charge contributions is 
ascertained. Sec. ^ deals with the spatial distribution of 
the WFs and the effect of molecular overlap on the di- 
electric properties of the compressed solid. In Sec. we 
investigate the nonuniqueness associated with the defini- 
tion of well-localized WFs in the dense solid. T he tig ht - 
binding analysis is presented in Sec. VII. In Sec. VIII we 
give a discussion of our results. 

Atomic units are used for all quantities except pres- 
sures, which are in gigapascal (GPa), and energies, in 
electron-volt (eV). Densities are expressed in terms of 
Ts, defined as (47rr^/3)aQ = V/N, where N is the num- 
ber of protons in the volume V, and oq is the Bohr ra- 
dius. Since the LDA tends to underestimate the pres- 
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where v is the volume of the unit cell of the crystal and 
the integral is over the Brillouin zone. Except for the con- 
straint V'n,k4-G = V'nk for all reciprocal lattice vectors G, 
the overall phases of the Bloch functions ■(/;„k = e*'' '"u„k 
are at our disposal. However, a different choice of phases 
(or "gauge"). 



Unk 



(2) 



does not translate into a simple change of the overall 
phases of the WFs; their shape and spatial extent will 
in general be affected, while the location of their centers 
of charge remains invariant modulo a lattice vector R 
If the band is isolated, Eq. ^ is the only allowed 
type of gauge transformation for changing the WF ?z;„(r) 
associated with that band. In the case of a composite 
group of bands, the allowed transformations are of the 
more general form 
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where Umn is a unitary matrix that mixes the bands at 
every wave vector k. Under this transformation the indi- 
vidual Wannier centers can shift, but their sum over the 
group of bands is preserved modulo a lattice vector 
Once a measure of localization has been chosen and the 
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group of bands specified, the search for the correspond- 
ing set of "maximally-localized" WFs becomes a problem 
of functional minimization in the space of the matrices 
The strategy of Ref. H is to minimize the sum of 



(k) 



the quadratic spreads of the Wannier probability distri- 
butions {|i(;„(r)|^}, given by 



(4) 



where the sum is over the chosen group of bands (in 
our application they will be the valence bands), and 
(r)^ = / r|w„(r)| dr, etc. Interestingly, the resulting 
"maximally- localized" (or "maxloc") WFs turn out to 
be real, apart from an arbitrary overall phase |3|. 

In numerical calculations the Bloch states -i/'nk are com- 
puted on a regular mesh of fc-points in the Brillouin zone; 
the integral in Eq. |l| is then replaced by a sum over the 
points in the mesh. In Ref. an expression was derived 
for the gradient of the spread functional H. with respect 
to an infinitesimal rotation SUml of the set of Bloch or- 
bitals, in terms of the Bloch functions in such a mesh. 
The only information needed for calculating the gradient 
are the overlaps (umk|u„.k+b), where b are vectors con- 
necting each mesh point to its near neighbors. Once the 
gradient is computed, the minimization can then proceed 
via a steepest-descent or conjugate-gradient algorithm. 

Since the Bloch eigenstates at different fc-points are 
initially computed by independent numerical matrix di- 
agonalizations, their phases are unrelated. As a conse- 
quence, the WFs obtained directly from them using the 
discretized version of Eq. |l] will be poorly localized, or 
not localized at all. In practice the following strategy 
is used for preparing a better set of Bloch states as the 
starting point for the minimization: one chooses a set of 
localized "trial functions" in the unit cell, which consti- 
tute a rough initial guess at the WFs; for solid hydrogen 
we use Gaussians on the centers of the molecules. Then, a 
unitary rotation among the initial Bloch orbitals is made 
in order to maximize their projections onto these trial 
functions (the detailed procedure is described in Eqs. 62- 
64 of Ref. 1^). For a reasonable choice of the width of 
the Gaussians (we have used a r.m.s. width of 1 A), the 
resulting rotated WFs are already extremely close to the 
"maxloc" ones, as discussed in Sees. VI and VII A| 



III. EQUILIBRIUM WANNIER DIPOLES 

The neutral entity composed of two paired nuclei and 
the occupied "maxloc" Wannier orbital centered around 
them forms a "Wannier molecule" in the bulk of the 
solid. In the low-density limit the "maxloc" WFs be- 
come nonoverlapping and coincide with the ground state 
bonding orbitals of isolated H2 molecules; however, at the 
high pressures we are interested in, there is an apprecia- 
ble overlap between neighboring WFs. In what follows 




FIG. 2. Contour plot of y'«Wi(r) for an occupied Wannier 
function in the Cmc2i structure at = 1.52 {v = 59.3 a.u. 
is the volume of the primitive cell). The central, cylindri- 
cally-shaped contour, which represents the bonding part of 
the WF, has a positive amplitude of +2.12; the outer lobes 
("orthogonality tails"), with antibonding character, have an 
amplitude of —0.11. 



we will sometimes loosely refer to the "maxloc" Wannier 
molecules in the dense solid simply as "molecules" . One 
should keep in mind, however, that had we chosen a dif- 
ferent measure of localization, the resulting "maximally- 
localized" WFs in the dense system would in general dif- 
fer somewhat from the ones we obtain. This nonunique- 
ness is intrinsic to WFs, and can be viewed as a manifes- 
tation of the ambiguity that always arises when trying to 
define "molecules" in a dense medium, either in terms of 
WFs or by oth er mean s. These issues will be discussed 
in Sees [yl and |VIIIB . 

Fig. 1^ shows a contour plot of a Wannier orbital for 
Cmc2i at = 1.52. The central positive contour with 
a large amplitude represents the molecular bond. The 
lowering of symmetry due to the crystalline environment 
is clear from the shape of the outer "corona" formed by 
the negative lobes, which have an antibonding character. 
These so-called "orthogonality tails" appear when the 
molecules overlap, due to the orthogonality requirement 
between different WFs; they are concentrated around the 
twelve nearest molecules, which allows for an efficient or- 
thogonalization between neighboring WFs. We will ar- 
gue in Sec. ^ that these overlapping orthogonality tails 
strongly influence the dielectric properties, and in partic- 
ular the vibron IR activity. 

An important effect of the anisotropic crystalline envi- 
ronment is that the molecules becomes polarized under 
the self-consistent internal electric fields inside the solid. 
In commonly used treatments of the dielectric proper- 
ties of molecular crystals |13,M, the so-called Clausius- 
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TABLE I. Static (Eq. 0) and vibron-induced (Eq. |) ef- 
fective charge vectors for the two molecules in the primitive 



cell of the 
symmetry. 


Cmc2i structure. 


The x-components vanish by 


Molecule 


Static 


In-phase Out-of-phase 


1 
2 




QyW qlW 9y(i) 

-4(1) 4(1) 4(1) -g°(l) 


TABLE IL Same as Table | 


|, but for the C2/m structure. 


Molecule 


Static 


In-phase Out-of-phaso 


1 
2 




4(1) ql{l) 4(1) g°(l) 

-4(1) -9^1) 4(1) 9°(i) 



Mossotti approximation is assumed: the system is mod- 
eled as a sum of nonoverlapping molecular charge distri- 
butions which become polarized in the local field pro- 
duced by the surrounding molecules; the bulk polar- 
ization is then the sum of the individual, nonoverlap- 
ping, molecular dipole moments, which can be straight- 
forwardly calculated from the bulk charge density p(r). 
Such a description becomes inappropriate whenever there 
is significant molecular overlap the electron den- 

sity becomes different from zero everywhere, and as a 
result the net dipole moment becomes dependent on the 
particular choice made for the unit cell [ p2p^ ]. In the 
case of molecular crystals, this is expected to occur, for 
instance, when the system is strongly compressed [ p4| ; 
under such circumstances a more careful treatment of 
the macroscopic polarization is required. According to 
the Berry phase theory of bulk polarization [ p3| , when- 
ever such overlap effects are significant, the macroscopic 
polarization Pmac of an insulating crystal cannot be ex- 
tracted from the bulk p{r), and is instead given, in the 
independent-electron approximation, by a Berry phase 
of the occupied Bloch states. This is a gauge-invariant 
quantity, and it is identical to another invariant, the to- 
tal sum of the dipoles d{n) of the "Wannier molecules" 
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V ^ 

n=l 



d(n). 



d(n) = -2e / r |w„(r)| dr, 



(5) 



(6) 



where M is the number of valence bands and e is the 
magnitude of the electron charge. In Eq. ^ the origin is 
chosen midway between the two paired protons, to cancel 
their contributions. The factor of 2 comes from spin- 
degeneracy (each occupied WF carries two electrons). 
We stress that Eq. || does not rely on the Clausius- 
Mossotti approximation at all, and it remains exact even 
when the "Wannier molecules" overlap strongly, as long 
as the system remains insulating. The decomposition 



of Pmac into individual Wannier dipoles is a powerful 
analysis tool, allowing us to go beyond the Berry phase 
approach used in Ref. |^ , which only gives the net po- 
larization of the unit cell. 

The permanent Wannier dipoles can be used to assign 
to each molecular charge distribution a "static effective 
charge" vector: 



q^(n) 



djn) 

rhond{n) ' 



(7) 



where rbond('T-) is the equilibrium bond length. This 
quantity, which vanishes in the low density limit of iso- 
lated molecules, measures the spontaneous symmetry- 
breaking charge transfer which occurs in the compressed 
solid whenever the two atoms in a molecule occupy crys- 
tallographically inequivalent sites. We emphasize that 
this definition is somewhat arbitrary, for the reasons dis- 
cussed at the beginning of this section (but see Sec. 0), 
and therefore (fin) does not relate directly to any mea- 
surable quantity. However, it is a sensible definition, 
which reduces to the natural one in the extreme ionic 
limit where the electron distribution is strongly concen- 
trated around the ions. In the next section we will use 
it to decompose the vibron effective charge into "static" 
and "dynamical" contributions, with the aim of under- 
standing the origin of the strong IR absorption in phase 
III. 

The location of the centers of the "maxloc" WFs re- 
flects the symmetry properties of the crystal; this is ap- 
parent from the form of the vectors q'^fri), shown in 
the first set of columns in Tables ^ and y for the two 
structures studied. The first set of columns in Tabic III 
lists their explicit values for = 1.52, both with the 
molecules on-site and after allowing them to move away 
from the ideal hep sites |l^j2^. For comparison we also 
report the values calculated using an electric quadrupolar 
(EQ) model |^^. In the C2/m structure there is a center 
of inversion between the two molecules in the primitive 
cell, leading to a cancellation of their permanent dipoles; 
in the lower symmetry Cmc2i the j/-components of the 
individual dipoles still average to zero over the primi- 
tive cell, but the z-components add up, yielding a small 
spontaneous polarization along the c axis (see Fig. |^). 
In general \qy \ » \ql\, so that the dipole moments make 
an angle with the molecular axes, i.e., in Table |l| is 
nonzero, although it is smaller than (7|. This agrees with 
the EQ model, where for the hep-centered (on-site) struc- 
tures the quadrupolar field at the center of the molecules 
is along y, so that the dipole moment along z is solely 
due to the small anisotropy in the polarizability |2^ . The 
EQ model also predicts larger dipoles in the C2/m than 
in the Cmc2i structure (although the effect is not nearly 
as pronounced as in the LDA WFs), as well as an in- 
crease in their magnitude as the molecules move off-site 
[|ll|,^. Although some of the qualitative features of the 
LDA Wannier dipoles are captured by the EQ model, its 
predictions are not reliable: for instance, it does not re- 
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TABLE III. Static (Eq. [^) and vibron-induced (Eq. |8|) effective charge vectors for molecule 1, for the hep-centered (on-site) 
and the off-site structures, at Vs — 1.52. Results are presented for the Wannier functions in the LDA approximation (WF) 
and for the electric quadrupole model (EQ). The i-components vanish by symmetry, gn and q± axe the magnitudes of the 
projections along the molecular axis and perpendicularly to it, respectively. 
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-0.003 
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-0.037 


0.161 


0.071 


-0.138 


-0.030 


0.129 


0.057 
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EQ 


-0.014 


-0.002 


0.012 


0.006 


-0.041 


-0.009 


0.038 


0.017 


0.004 


-0.002 


0.002 


0.004 
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0.019 
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-0.043 
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0.004 
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0.073 


0.013 
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-0.437 
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0.058 


0.013 
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WF 
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-0.054 
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0.115 
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FIG. 3. Molecular effective charge vectors for molecule 
1 in the hep-centered (on-site) Cmc2i structure at several 
densities. • Static charge (Eq. R); B dynamical component of 



effective charge (last term in Eq. ^ for the in-phase vibron; ▲ 
dynamical component of effective charge for the out-of-phase 
vibron. The i-components vanish by symmetry. 
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FIG. 4. Same as Fig. q, but for the hep-centered (on-site) 
C2/m structure. We have only considered densities for which 
the LDA band gap remains open, which are below the density 
at which phase III appears. 



produce the change in sign of ql{n) for Cmc2i at low 
pressures other discrepancies can be seen in Table 
[II, most notably in th e vibron effective charges, and are 
discussed in Sec. |IVB . 

From Fig. ^ we can already see that |q''(n)| will be 
small, since the Wannier distribution is fairly symmetric 
with respect to the center of the paired protons. The 
dependence of the static charges versus is plotted in 
Figs. ^ and 0; as expected they vanish in the low den- 
sity (large Vg) limit, and even at the highest pressures 
(~ 210 GPa) they are only a few percent of the electron 
charge, indicating that, at least in the structures under 
consideration, the ionicity of the molecules remains quite 
small, contrary to some proposals |27|j2^ ]. Notice also 
that at high pressures (f{n) becomes quite sensitive to 
the crystal structure (see Table HI). This is not surpris- 



ing, since it is totally induced by the crystal field. At 
^ 165 GPa {vs = 1.45) the permanent dipole moment 
of an H2 Wannier molecule in the hep-centered Cmc2i 
structure becomes 0.075 a.u., i.e. more than 1/10 of the 
dipole of an isolated water molecule (0.74 a.u.). 



IV. VIBRON INFRARED ACTIVITY 

A. Wannier-function description 

The investigation of the number of IR-active lattice 
modes and their oscillator strengths in candidate struc- 
tures is a useful guide in the search for the structures of 
the high-pressure phases |l^,|l^,^ . Here we will focus on 
the vibron IR absorption observed in phase III llqjl^] • 
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relatively strong intensity is somewhat puzzling, since the 
stretching mode of the isolated H2 molecule is Raman- 
active but IR-forbidden. and this has stimulated a large 
number of studies 
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IR absorption is caused by the coupling of light to 
the change in Pmac induced by the lattice modes. In 
the basic theory of IR absorption in molecular crystals 
|^ , pT| , the Clausius-Mossotti approximation of nonover- 
lapping molecules is assumed. The modern theory of 
polarization |p3| treats rigorously the situation where 
that approximation breaks down, as discussed in the pre- 
vious section. There, we decomposed the spontaneous 
macroscopic polarization into contributions from individ- 
ual Wannier molecules; here we will do the same for the 
vibron-induced fluctuations in Pmac. For every vibron 
mode v we will assign to each Wannier molecule a vibron- 
induced "effective charge" vector (compare with Eq. |^): 



dd{n) 



(8) 



where Ui, is the normal coordinate [ p^ . The above ex- 
pression is very similar to the definition of the Born ef- 
fective charge associated with the stretching mode of a 
diatomic molecule |Q. The vector q'^(n) measures the 
vibron-induced symmetry-breaking charge transfer; with 
the help of Eq. it can be decomposed into two parts 
(see Eq. 20 of Ref. |33|) pi: 



d(f{n) 
dui, 



(9) 



Since (f{n) is the static charge, we will call the second 
term on the r.h.s. the dynamical charge. 

Experimentally two vibrons have been detected in 
phase III: the lower frequency mode appears in the Ra- 
man spectrum, and the higher-frequency one in the IR 
spectrum [p|jl5[. Both the C2/m and the Cmc2i struc- 
tures have two vibron modes: one in which the two 
molecules in the primitive cell vibrate in-phase (i^ = i), 
and a higher frequency mode in which they vibrate out- 
of-phase (z/ = o) . The effective charges q"^ (n) were calcu- 
lated using Eq. ^ by changing the molecular bond length 
by small amounts 6ui, in the range [0.0015,0.0035] a.u., 
after checking that such displacements yield essentially 
linear changes in the Wannier dipoles. Table III shows 
their values for Vg — 1.52; as in the case of the static 
charges, we have in general that gj^ > g^. Figs. ^ and 
^ plot the static and dynamical charges versus . Since 
these originate from the interactions between molecules, 
their magnitudes vanish in the low-pressure (large r^) 
limit, and increase as pressure goes up. The most striking 
feature is that the dynamical terms increase with pressure 
much more rapidly than the static ones; at the highest 
pressures studied they are already 3.3 to 6.7 times larger, 
depending on the structure and on the vibron mode. 
Since this appears to be a rather general feature, it is 
also likely to occur in the yet-undetermined structure of 



phase III; the observed strong vibron IR activity is prob- 
ably caused by this increase of the dynamical charges. 
Their dominant role had been previously inferred from 
the strong anisotropy of the atomic Born effective charge 
tensors [E5|. 

Figs. I^and ^ and Table III also show that the vibron- 
induced fluctuations in the individual molecular dipoles 
(q'^(n)), although clearly mode-dependent, are compa- 
rable for the two vibrons. (I nteres tingly, this is not so 
for the EQ model-see Section IV B.) The important dif- 
ference occurs only after adding the contributions from 
the two molecules, and can be seen in Tables | and ||. 
In the in-phase mode the large y-components cancel be- 
tween the two molecules in the primitive cell, resulting 
in a weak IR activity (which actually vanishes in C2/m, 
since the small z-components also cancel). By contrast, 
in the out-of-phase mode the large y-components add up, 
resulting in a large net dP^^c/ du^, and thus in a strong 
IR activity. 

We emphasize again that all the quantities in Eq. || are 
gauge-dependent, like the Wannier functions themselves 
(but see Sec. VI). The gauge-invariant, measurable quan- 



tity is the net "vibron effective charge" vector, obtained 
by averaging q"^ (n) over all molecules in a primitive cell: 



M 



(10) 



where nmoi is the number of molecules per unit vol- 
ume; the vibron oscillator strength is proportional to 
'T-moilq'^l^- The calculated values of |q''| versus are 
plotted in Fig. ^ together with the experimental results. 
For the out-of-phase mode in the Cmc2i structure the 
LDA calculation yields values very close to the exper- 
imentally measured IR absorption in phase III, but on 
the other hand the IR activity of the in-phase mode, 
although weaker, would still be observed, which is not 
the case. As for C2/m, the IR absorption is too strong 
compared to experiment. Hence it seems that neither 
structure is likely to be the correct one for phase III (this 
is also supported by the large number of observed libron 
modes in phase III p5| , ^ , which is incompatible with 
structures with such small primitive cells). Nevertheless, 
the above results for these structures allow us to make 
an important general point. They show that large per- 
manent molecular dipoles are not required in order for 
strong vibron IR absorption to occur, contraryto what 
has been sometimes stated in the literature j|l^,^,0. 
In fact, in both structures the magnitude of the perma- 
nent dipoles (static charge) is far too small to account 
by itself for the measured absorption. However, once the 
dynamical charge transfer is accounted for, the resulting 
IR activity becomes even larger than the one measured 
in phase III. 

Finally, Table III and Fig. |^ show that the displace- 
ment away from the hep sites significantly increases the 
vibron-induced charges, as well as the static charges. 
This is yet another example of the strong sensitivity of 
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■ Cmc2| (in-phase) 

i Cmc2j (out-of-phase) 
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FIG. 5. Magnitude of the net vibron-induced effective 
charge vector Iq"! (Eq. ^) versus r^, for the in-phase and 
out-of-phase vibron modes in the Cmc2i and C2/m struc- 
tures. Experimental data for phases II and III is from 
Ref. 11 7|, and was converted from Szigeti to Born charges |p8|. 



the effective charges to the crystal structure, which may 
help explain the large difference in the intensity of IR ab- 
sorption between phase III and the lower pressure phases 

p5|,ra. 



B. Comparison with the EQ model 



Table [11 also contains the values of the vibron-induced 



charge vectors on each molecule, as calculated from the 
EQ model |2^. Like the static charges, they are signif- 
icantly smaller than those obtained from the WFs. An- 
other important difference is that in the EQ model the 
vibron-induced charges on the individual molecules are 
more than one order of magnitude smaller for the out- 
of-phase than for the in-phase mode, whereas their WF 
counterparts are comparable for the two vibrons. The 
reason is the following: in the EQ model the vibron- 
induced change in a molecular dipole can be written to 
first order as 



6d\\ c± 5a\\E\\ + ol\\SE\\ 



(11) 



where a is the molecular polarizability and E is the 
quadrupolar electric field on the molecular site (for defi- 
niteness we look at the dipole along the molecular axis; 
the same analysis applies to the perpendicular compo- 
nent). The first term on the r.h.s. is equal for the 
two vibrons, so that the difference between their effec- 
tive charges arises from the second term. Choosing the 
isolated-molecule parameters from Refs. |^ and it 
turns out that the two terms have a very similar magni- 
tude. But whereas in the in-phase mode they have the 
same sign, in the out-of-phase mode they have opposite 
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FIG. 6. Modulus squared (upper panel) and dipole (lower 
panel) of a WF, accumulated by integrating up to a certain 
radius around the molecular center, for Cmc2\ at = 1.52. 
The a;-component of the dipole vanishes by symmetry inde- 
pendent of radius, and the horizontal arrows denote the con- 
verged values. The thin lines correspond to the WF obtained 
from bond-centered Gaussians with a r.m.s. width of 2.0 A 
(see Sec. |l]). 

signs, so that their contributions largely cancel, resulting 
in a much smaller molecular effective charge. As a con- 
sequence, in the Cmc2i structure the in-phase oscillator 
strength comes out larger than the out-of-phase, which 
is the opposite of the LDA result. These discrepancies 
between the LDA WFs and the EQ model are likely to 
be related at least in part to the rather delocalized na- 
ture of the induced dipoles, which will be discussed in the 
next section, whereas the EQ model assumes point-like 
(infinitely localized) molecules. 



V. SPATIAL EXTENT OF THE WANNIER 
MOLECULES 

In the previous sections we focused our attention on 
information that can be extracted from the location of 
the centers of charge of the WFs. Here we will exam- 
ine in detail their spatial distribution at high pressure, in 
particular their spatial extent. This will allow us to in- 
vestigate the effects associated with the overlap between 
neighboring WFs; such effects are expected to be sig- 
nificant at megabar pressures, as suggested by the large 
bandwidths, in excess of 20 eV (see Fig. n3). 
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FIG. 7. Fraction of the Wannier dipole dy{l) accumulated 
by integrating up to a certain radius around the molecular 
center, versus the fraction of the Wannier charge that lies 
inside the same radius, for Cmc2i (see also Fig. H). 



A. Spread of the Wannier charge and dipole 
distributions 

The spread of the Wannier charge and dipole distribu- 
tions are presented in Fig. ^. For Cmc2i at = 1.52 
(r^ = 2.0), a radius of 1.72 a.u. (2.25 a.u.), half the 
shortest intermolecular distance (iH2-H2j encloses about 
72% (85%) of the charge and only 17% (28%) of the y- 
component of the dipole. This suggests that already at 
Vs = 2.0 (~ 13 GPa) the overlap between nearby WFs is 
far from negligible. It is also clear that the dipole is signif- 
icantly more spread out than the charge, with very large 
contributions arising from the orthogonality tails in the 
overlap region, where the Wannier charge density is very 
small. The longer range of the dipole is to be expected, 
due to the factor r in Eq. but the large magnitude 
of the effect is somewhat surprising. Fig. |^ shows even 
more clearly that for both pressures a rather small frac- 
tion of the total charge, located in the overlapping tails 
of the WF, is responsible for most of the dipole. Also 
striking is the fact that, at = 1.52, up to a radius of 
7 a.u. the accumulated (iz(l) remains positive, whereas 
the converged value is negative; this suggests that the 
agreement in sign with the EQ model (see Table III) may 
be fortuitous, since in that model the dipole is caused by 
the electric field at the center of the molecule. 

By analogy with the radially integrated Wannier 
charge and dipole distributions (Fig. one can plot the 
derivative of these quantities with respect to the normal 
coordinate of a vibron mode (Fig. ^). At high pressures 
the contributions from the overlapping tails to the change 
in the dipole moment are very significant, even more so 
than for the equilibrium dipole. In other words, in the 
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FIG. 8. Upper panel: derivative with respect to the 
out-of-phase vibron's normal coordinate Uo of the accumu- 
lated probability plotted in the upper panel of Fig. |^. Lower 
panel; derivative with respect to Uo of the accumulated dipole 
plotted in the lower panel of Fig. |^, i.e., accumulated radial 
integral of the vibron-induced effective charge vector q°(l) 
(Eq. 1^ (the a;-component vanishes by symmetry for all ra- 
dius). The arrows denote the converged values. 



dense solid the dynamical charge transfer processes re- 
sponsible for the IR activity are very delocalized. 

The effect of a vibron on the Wannier charge distri- 
bution is depicted in the upper panel of Fig. 0: charge 
is depleted from the inner part and accumulates in the 
outer part of the molecule. Note that the charge trans- 
fer occurs mainly along the molecular axis, and is essen- 
tially symmetrical with respect to the molecular center, 
as one would expect from stretching an isolated molecule, 
which does not break the symmetry between the two 
atoms. This kind of charge transfer alone would lead 
to a zero net change in the (vanishing) molecular dipole, 
and hence to no IR absorption. The contribution of a 
Wannier molecule to the IR activity of the crystal comes 
from the comparatively small odd part (with respect to 
its center) of the charge transfer. Since this is barely 
visible in the upper panel of Fig. |^, in the lower panel 
we have removed the large even part. It is interesting to 
note that near the paired protons most of the odd part 
is oriented roughly perpendicularly to the molecular axis, 
making a small angle with the c axis, such that it gives a 
small positive contribution to ^^^(l) and a larger negative 
contribution to q°{l). This observation is supported by 
Fig. |8[ which shows that for small radius the accumulated 
q°{l) is negative and larger than the accumulated 9^(1), 
which is positive. For large radius 9^(1) changes sign and 
ends up overtaking 9°(1), and the net molecular vibron 
charge vector has a larger projection along the molecular 
axis than perpendicularly to it (see Table III). 
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The results of this section should be relevant for mod- 
els that attempt to account for the dielectric properties 
of compressed hydrogen. For instance, it seems unlikely 
that models based on point-like objects, such as the elec- 

contain all the 
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trie quadrupole (EQ) model 
important ingredients that lead to the strong IR absorp- 
tion in the highly compressed phase III. The reasons are 
twofold: at such high densities (i) the electrostatic inter- 
actions are expected to differ substantially from the ideal 
quadrupolar one ||4l|] (and in fact the validity of a multi- 
pole expansion becomes questionable when the molecu- 
lar charges overlap significantly), and (ii) the "classical" 
treatment of polarization, based on the bulk p{r), be- 
comes inadequate [|3|. We note that although the EQ 
model can account f or bot h static and dynamical charges 
||l|,||j3l) (see Sec. |VB|), both effects are then due to 
local fields and polarizabilities (see also Fig. 2 of Ref. js^ j 
and associated discussion regarding local versus nonlocal 
mechanisms) . 

The central conclusion of the preceding analysis is that 
the contributions to the induced molecular dipoles (and 
their fluctuations) arising from the overlapping tails of 
the Wannier orbitals, which extend well beyond the near- 
est neighbor molecules, are crucial. It is instructive to 
contrast this state of affairs with what happens in liq- 
uid water: there, the contribution from the orthogonal- 
ity tails is negligible . This difference may stem from 
the fact that an isolated water molecule is polar, so that 
the effect of the liquid environment is only to modify a 
previously existent dipole moment, whereas in solid hy- 
drogen the molecular dipole is totally induced. Induced 
dipoles tend to be rather extended because the outer re- 
gions of the molecules are most easily polarizable [ p3[ ; 
thus their contribution to the dipole can be large, even 
though p(r) is small, because of the r factor in Eq. |^. 
In conclusion, the WF analysis strongly suggests that 
the Clausius-Mossotti picture of nonoverlapping dipoles 
breaks down rather dramatically for solid hydrogen at 
megabar pressures. 



B. Measuring the molecular overlap 



FIG. 9. Upper panel: derivative of the modulus squared 
of the Wannier orbital with respect to the normal coordinate 
of the out-of-phase vibron mode {vd\wi(r)\^ /duo) for Cmc2\ 
at rs — 1.52 {v — 59.3 a.u. is the volume of the primitive 
cell). The central contour has an amplitude of —4.5, and the 
two outer contours have an amplitude of -1-0.25. Lower panel: 
odd part (with respect to the center of the molecule) of the 
same quantity. The upper (lower) contour has an amplitude 
of +0.05 (-0.05). 



The overlap between the charge distributions of neigh- 
boring WFs can be quantified as @J42) 

^ J\wm{v)\^\wn{r)fdr 

^mn — -^/2 1/2' V / 

[jMr)\'dr) (/|«;„(r)|*dr) 

For Cmc2i the largest value of Omn is 0.005 at — 2.0 
and 0.021 at rg = 1.52; the latter value is still quite small, 
roughly twice the value for WFs located on nearby water 
molecules in liquid water Thus, by inspection of 

Omn alone one would not suspect that the overlapping 
tails are so much more important for the dipole moments 
in compressed solid hydrogen than in liquid water. The 
reason is that Omn measures the overlap between charge 
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FIG. 10. Static charge vector q''(l) (Eq. 0) versus n, where 
n X n X n is the size of the mesh of fc-points used for computing 
the WFs, for the on-site Cmc2i at Va = 1.52. The number of 
fc-points is kept fixed at 11x11x11 during the self-consistent 
calculation. 



distributions, whereas in this system the dipoles are much 
more spread out. 

Another indication that overlap effects are important 
comes from the well-known fact that a very large num- 
ber of fc-points is required to conver ge t he total-energy 
calculations in compressed hydrogen [ [ll| . In fact, if the 
molecules were strictly nonoverlapping a single /c-point 
would suffice for computing all physical properties. In 
Fig. |l^ is shown the static charge (i.e., the dipole) for 
WFs obtained using different meshes of A:-points. It is 
clear that a dense mesh is required for converging this 
quantity. This results from the fact that, when using 
a discrete mesh, the WFs are actually periodic in real 
space, with a periodicity which is inversely proportional 
to the spacing between neighboring points |^ . Hence the 
need for a fine sampling of the Brillouin zone is just a 
manifestation of the large contributions to the Wannier 
dipole arising from the tails far away from the "home" 
unit cell. In solid hydrogen a dense mesh of /c-points 
is expected to be even more important for the dielectric 
properties than for the total energy, since Pmac is partic- 
ularly sensitive to the Wannier tails. 
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FIG. 11. Root-mean-square width A = {D,/3M) ' of the 
WFs and electronic localization length ^, both averaged over 
all directions. The dashed line denotes half the shortest in- 
termolecular distance in Cmc2i (in C2/m it is very similar). 



els of molecular solids [^Q. This is an overlap effect, 
caused by the orthogonality requirement. It is due to 
the enhancement of the outer corona shown in Fig. ||, 
which is not included in the definition of "molecules" in 
those models. It can also be viewed as a result of the gap 
reduction with pressure. 



Also plotted in Fig. [11 
length I = (f7i/3Af)^/^ 



is the electronic localization 
[lHJill, where Qi < 



Q is the 

gauge-invariant part of the spread of the WFs ^ 
measures the mean-square quantum fluctuation of the 
macroscopic polarization [Q, normalized in such a way 
as to be finite for insulators, and diverging when the band 
gap closes. Notice that in the low-density limit A — > ^; 
that happens because there is only one occupied WF per 
molecule, and can be understood by comparing Eqs. 14 
and 15 of Ref. As expected, at high pressures A and ^ 
are larger for the C2/m structure, which has the smaller 
band gap: at — 1.52 A has increased by 13% (20%) in 
Cmc2i (C2/m) with respect to the low-density (isolated 
molecule) value. 



C. Quadratic spread and localization length 

Another way of quantifying the spatial extent of the 
WFs is in terms of their quadratic spread. According to 
Eq. ^, the r.m.s. width of the Wannier probability distri- 
bution, averaged over all three Cartesian directions and 

— 1/2 
over the occupied WFs, is A = {Q/3M) . This quantity 

is plotted versus in Fig. |l]. Notice that it increases 
with increasing pressure (decreasing Tg), i.e., the Wan- 
nier molecules become more extended upon compression, 
which is the opposite of what happens in the usual mod- 



VI. UNIQUENESS OF WELL-LOCALIZED 
WANNIER MOLECULES 

As discussed earlier, the WFs are strongly nonunique; 
in particular, it is only the sum of all the Wannier 
dipoles over a primitive cell which is physically mean- 
ingful (Eq. 1^) , whereas the individual dipoles (and hence 
cf{n) and q^{n)) are gauge-dependent. Nevertheless, in 
this work we have been looking at the individual dipoles 
of the "maxloc" WFs in an attempt to extract from them 
useful physical information. The underlying assumption 
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is that in practice well-locahzed WFs are fairly unique. 
For solid hydrogen this is obviously true in the low den- 
sity limit, where they reduce to the bonding orbitals of 
the individual molecules. Here we will discuss to what 
extent that assumption holds for the compressed solid as 
well. 

A systematic way of assessing the degree of unique- 
ness of well-localized WFs would be to implement dif- 
ferent localization criteria and then compare the result- 
ing "maximally-localized" WFs. We have not attempted 
such a detailed study; instead we have performed a sim- 
pler test, which is a first step in that direction. As men- 
tioned in Sec. in the method we are using |^ an ini- 
tial guess is made for the localized WFs, with the help 
of "trial functions" , which in our case are bond-centered 
Gaussians (we will call the resulting orbitals "projected 
WFs" ) . Their localization is then enhanced by minimiz- 
ing the quadratic spread fl (Eq. ^ , yielding the "maxloc" 
WFs. Since the projected WFs are totally oblivious of 
the localization criterion that one later uses to further 
localize them, it seems reasonable to assume that the 
difference between the projected and the "maxloc" WFs 
is an upper bound to the differences that would occur 
between WFs obtained using any two "sensible" local- 
ization criteria. 

Let us consider the Cmc2i structure at = 1.52, for 
which the average r.m.s. width of the "maxloc" WFs is 
A = 1.11 a.u. If we choose the r.m.s. width of the initial 
bond-centered Gaussians to be 1.89 a.u. (1 A), the result- 
ing projected WFs are essentially indistinguishable from 
the "maxloc" ones: for instance, the curves correspond- 
ing to those in Fig. ^ are virtually identical, and the in- 
dividual Wannier dipoles remain the same to at least six 
significant digits! This is compelling evidence for a high 
degree of uniqueness of well-localized WFs in this sys- 
tem, at least for the high-symmetry configurations that 
we studied. If we double the width of the initial Gaussian, 
some differences start to appear. They remain barely vis- 
ible in the accumulated radial integral of the probability 
(upper panel of Fig. ||) but are noticeable, although still 
relatively small, in the radially integrated dipole (lower 
panel of Fig. For instance, the large y-component of 
the dipole changes by around 2%. 

VII. TIGHT-BINDING ANALYSIS 

In this section, we investigate whether the essential 
physics of the WFs in high-pressure H2 phases can be 
captured by a simpler tight-binding (TB) approach. As 
is well known, the TB approximation provides a simple, 
computationally inexpensive method for computing elec- 
tronic structure effects, and has the additional advantage 
that its output is easily interpreted in terms of a local, 
real-space picture . Thus, TB is a natural approach to 
explore here, where we want to study the dielectric prop- 
erties of H2 phases from just this kind of local point of 



view. We confirm below that a previously-proposed sp^ 
TB model |Q provides a good description of the occu- 
pied bands in these systems, and show how the operations 
of constructing the WFs and computing their contribu- 
tions to dielectric properties (such as electric polariza- 
tion) can be carried out in the TB framework. Finally, 
using the fact that the TB representation automatically 
provides an atom-by-atom and orbital-by-orbital decom- 
position, we obtain useful insights into the nature of the 
WFs and their contributions to the dielectric properties. 

A. Tight-binding formalism 

In the TB method, the Bloch functions Vnk are ex- 
panded in a basis of atomic-like orbitals (f)ii as 

V'nk(r) =^C„k(«Oe'*'""'/'.i(r) ■ (13) 

il 

Here I labels the unit cell located at R/, i labels an or- 
bital on the atom at r^; = Ri + Ti (where specifies the 
relative position of the atom within the unit cell), and 
the vector of coefficients C„k(iO forms the TB represen- 
tation of the Bloch function. Our goal is to carry out a 
unitary transformation to a set of M localized WFs 

= ^ E^-^"k(r) (14) 

nk 

associated with a set of M occupied bands, where the 

(k) 

Uan are k-dependent M x M unitary matrices that will 
be fixed by the requirement of maximal localization . 
Introducing the TB representation of the WF 

Wa{r) =J2w<^{il)My^) > (15) 

il 

it follows that 

Wo.m = ^ E e''"'"" ^»k(*0 . (16) 

nk 

The essential ingredients needed for the construction of 
the maximally-localized WFs , or for the computation 
of the Berry-phase polarization ||2^, are inner products 
(w„k|u„'k') between the cell-periodic part of the Bloch 
functions 

wnk(r) = e-'''"-^„k(r) (17) 

at nearby fc-points in the Brillouin zone. In principle, the 
calculation of the (wnk|un'k') requires a detailed knowl- 
edge of the basis orbitals (pu, which has been done in 
Ref. |]5l|] . However, in the spirit of minimal empirical 
TB, we make the approximation 

{Un^lUn'i.') = E^"k(*0 C„,k'(*0 ■ (18) 
il 
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[When completing the circuit across a Brillouin zone 
boundary, the relation C„^k+G(*0 = e^'*^ '''' C„k(iO 
should be used to translate by a reciprocal lattice vec- 
tor G.] 

Eq. can be derived via a Taylor expansion of the 
exponential factor exp[ik • (r — r^;)], with the following 
assumptions: (i) that the TB basis orbitals are orthonor- 
mal, {(j)ii\4>i'i') = Sii'dw; (ii) that the position operator 
is diagonal in the TB basis, {(l)ii\'c\4>i'i') — vu 6^ Sw; and 
(iii) that matrix elements of higher powers of the position 
operator are likewise trivial. 



TABLE IV. Tight-binding parameters of Ref. pd, in eV. 



■u\xPy''2 



q r 
Vil Z^l 



(19) 



Conditions (i) and (ii) are actually special cases of (iii) 
and all all are quite artificial in that they cannot be satis- 
fied for actual basis functions. For example, while an sp^ 
hybrid on a given atom should have its charge center dis- 
placed from the geometric center of the atom, condition 
(ii) does not allow this effect to be captured. Similarly, 
the spread {<t'ii\r'^\4>ii) ~ {4'ii\'^\4'ii)'^ of individual ba- 
sis orbital is taken to vanish, according to condition (iii) . 
Nevertheless, Eq. (|l^) is the logical extension of the em- 
pirical TB philosophy, in which one tries to avoid intro- 
ducing any additional parameters beyond those needed 
to parameterize the Hamiltonian itself. Despite its sim- 
plicity, the tests presented below demonstrate that this 
approach captures much of the interesting complexity of 
the WFs, at least for the systems under study here. A 
similar approximation was previously shown to allow for 
reasonably accurate TB calculations of dynamical effec- 
tive charges in semiconductors . 

In practice, we work entirely within the TB representa- 
tion. First the C„k(*0 ^-re determined on a regular mesh 
of fc-points by solving the standard secular equation in- 

fk) 

volving the Hamiltonian matrix H^, . Then the electric 
polarization can be cornputed by inserting Eq. ( |l8| ) into 
the formalism of Ref. |23|] . Similarly, an "optimal" set 

(k) 

of unitary matrices Uan can be obtained by inserting 
Eq. (|8|) into the formalism of Ref. 10, and from these, 
the WFs Wa{il) obtained via Eq. (p^. The resulting 
WFs are optimal in the sense of being maximally local- 
ized in real space, i.e., of minimizing Eq. (1). 

In the "maxloc" method one usually begins by 
choosing a set of localized "trial functions," and mak- 
ing a preliminary unitary rotation among the Bloch or- 
bitals in order to maximize their projections onto these 
trial functions, as discussed in the LDA context at the 
end of Sec. ||. In the TB context, we have found the 
following natural way of constructing the trial Wannier 
functions. Since we have two molecules per cell, we want 
to carry out the 2x2 rotation that makes one state have 
most of its projection on the first molecule, and the other 
have most of its projection on the second molecule. To 
do this, we consider the difference AP = Pi — P2 of 
projection operators Pi and P2 onto the first and sec- 
ond molecule, respectively. (In the TB basis, AP is just 
a diagonal matrix with ±1 diagonal entries.) Then, at 



Hopping parameters 


Intramolecular 


Intermolecular 


Vss 


-8.50 


-0.04 




-8.75 


-0.16 




+9.00 


+0.89 



each k, we diagonalize AP in the space of the two Bloch 
states, and set the phase of each eigenvector by requiring 
that its inner product with an even linear combination 
of s orbitals on the two atoms comprising the molecule 
should be real and positive. We find that the unitary 
transformation Uan obtained in this way turns out to be 
an excellent approximation to the "maxloc" Um which 
minimizes Eq. (1). In fact, subsequent minimization typ- 
ically only leads to changes of WF coefficients of order 
one part in 10~^, and so that in practice it is not even 
necessary to carry out the maxloc minimization proce- 
dure. This is consistent with our similar experience in 
the LDA context as discussed at the end of Sec. Vl. 



B. Details of the tight-binding model 

The TB parameterization we used is the one proposed 
by Chacham et al. [Q. These authors showed that the 
main characteristics of the electronic structure of high- 
density solid hydrogen at megabar pressures could be re- 
produced by using a minimal orthogonal TB basis com- 
prised of s, Px, Py, and p^ orbitals on each hydrogen 
atom. The intermolecular matrix elements are taken as 

= K.,(do)e"(i-''/*) 

Kp(do)et("+'')/'l(i~'^/*) 
(20) 



rpp, 
Vsp 



with dimensionless constants a = 5.76 and /3 = 2.52, 
where d is the interatomic distance and c?o — 3.79 A is the 
equilibrium hep lattice constant at zero pressure. Given 



in Table IV are the Vss (do), Vppa-{do), Vsp{do), and the 



intramolecular matrix elements Vss, Vs, 



and Vppa which 



are independent of d. The intra-atomic parameter Cs — 
ep = -20 eV. 

In addition to this original tight-binding scheme, we 
have also tested an extended scheme that includes a 
correction designed to incorporate the effects of the 
quadrupolar electrostatic fields arising from neighbor- 
ing molecules. In this extended "TB-I-Q" scheme, the 
molecules are first modeled as point quadrupoles cen- 
tered at the molecular sites (mid-bond positions). The 
quadrupole moment tensor for each molecule is taken 
from the free-molecule calculations of Ref. |^9| by as- 
suming a linear dependence upon the bond length in the 
range of 1.4-1.6 a.u. The total quadrupolar electric field 
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is then evaluated at each molecular site, and the electro- 
static potential shift on each atom in the molecule is cal- 
culated by assuming a linear extrapolation to the atomic 
position. Finally, the diagonal elements (self-energies) 
of the TB Hamiltonian matrix are modified by adding 
these energy shifts, and the solution of the secular equa- 
tion then proceeds as usual. 



C. Band structure 

We have applied the TB model to the same Cmc2i 
and C2/m candidate H2 structures studied with LDA 
methods in earlier sections, and confirmed that this TB 
model does a good job of reproducing the critical fea- 
tures of the electronic band structure. For convenience, 
we present only results on the Cmc2i geometry; the cor- 
responding results for C2/ m a re qualitatively similar. As 
already indicated in Sec. IC, our Cmc2i structure has 
Ts — 1.52a.u., Tbond = 1.445a.u, c/a — 1.576, and the 
tilt angle 6 — 54.0°. The two-molecule (four-atom) unit 
cell is illustrated in Fig. |l|. Use of the sp'^ TB basis leads 
to a 16 X 16 TB Hamiltonian matrix. 

Figure |l^ shows the good agreement between TB and 
LDA band structures for this geometry. The agreement 
in the occupied valence-band region (lowest two bands) 
is excellent, and the resemblance in the conduction-band 
region is also reasonable. The TB model predicts a gap 
closure at a density of 0.3962 mol/cm"^, consistent with 
the results from other studies |5^. The band structure 
is hardly affected at all if the TB+Q theory is used in 
place of simple TB. 

The WFs are constructed for the two occupied bands 
by using the definition in Eq. |l^. A 10 x 10 x 10 /c-point 
mesh is used in calculation. Since there are only two 
WFs in the unit cell, and these are related to each other 
by a symmetry, it suffices to analyze just one of them. 
In Table ^ we analyze the spatial distribution of the 
WF by decomposing into contributions coming from the 
"home molecule" and the first two nearest-neighbor shells 
of molecules in real space. The home molecule is labeled 
as "Neighbor 0" , the next six neighboring molecules form 
the first shell at a radius of 3.4342 a.u. (0.9767 in units of 
lattice constant), and the second shell comprises the next 
six neighboring molecules at a radius of 3.5162 a.u. (one 
lattice constant) . Table ^ gives a clear picture of the spa- 
tial structure of the WFs. We can see that the WFs have 
about 85% of their probability on the home molecule, 
91% (cumulatively) inside the first shell, and 97% up to 
the second shell. In the TB framework, the contribution 
to the WFs can be very easily decomposed further into 
the sp3 tight-binding basis orbitals, as shown in Fig. 
While the s-orbital contribution is ~5 times larger than 
that of the p orbitals, we find that the s contribution is 
almost entirely localized to the home molecule. On the 
other hand, although the p orbitals give a smaller total 
contribution, they play a much bigger role in the tail re- 
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FIG. 12. Electronic band structure calculated by (a) 
tight-binding method, and (b) LDA approach. 



gion that determines the spatial distribution of the WFs. 



D. Tight-binding Wannier functions 

The above quantities have been recomputed using the 
TB+Q theory in place of the simple TB theory, but the 
differences are not significant. 

According to the modern understanding , the elec- 
tronic contribution to the polarization can be equiva- 
lently expressed either in terms of a Berry phase (BP) 
computed from the Bloch functions, or in terms of the 
displacements of the Wannier centers (i.e., from the 
molecular dipole moments). These two approaches are 
compared for each of the three different computational 
schemes (TB, TB+Q, and LDA) in Table 0. It can be 
seen that the bulk polarizations Pz computed from the 
WF and BP approaches are in good agreement with each 
other (^5%) for all three schemes, indicating good inter- 
nal consistency. (The small discrepancies can be traced 
mainly to incomplete fc-point convergence.) 

The advantage of the WF approach is that the decom- 
position into molecular dipole moments can give some in- 
sight into the microscopic origins of the dielectric prop- 
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TABLE V. Spatial distribution of Wannier functions for 
the first three shells. Shown in the table are the list of the 
nearest neighboring molecules, their distance to home unit 
cell, the contribution to probability from each molecule, and 
the accumulated probability up to the current molecule. 



Neighbor 


Radius (a.u.) 


Probability 


Accum. Prob. 





0.0000 


0.84337 


0.84337 


1 


3.4342 


0.00764 


0.85101 


2 


3.4342 


0.00764 


0.85865 


3 


3.4342 


0.01294 


0.87159 


4 


3.4342 


0.01294 


0.88453 


5 


3.4342 


0.01898 


0.90352 


6 


3.4342 


0.00648 


0.91000 


7 


3.5162 


0.00984 


0.91984 


8 


3.5162 


0.00984 


0.92969 


9 


3.5162 


0.01098 


0.94067 


10 


3.5162 


0.01098 


0.95165 


11 


3.5162 


0.00777 


0.95942 


12 


3.5162 


0.00777 


0.96719 
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FIG. 13. The spatial distribution and decomposition of 
Wannier functions. O is the contribution to probability from 
p-orbitals, O is the contribution from the s-orbitals, and (A) 
is the sum of the s and p contributions. 



TABLE VI. Comparison of dipole moments and bulk po- 
larization calculated by TB, TB-I-Q and LDA. Because of the 
symmetry of Cmc2i, dix = d2x = 0, diy — —d2y (so Py — 0), 
and diz — d2z- The results from the Wannier- function analy- 
sis are given in the first three rows, while the bulk polarization 
from Berry-phase calculation appears in the last row. Atomic 
units are used. 





TB 


TB+Q 


LDA 




-0.03326 


-0.03759 


-0.06297 


diz 


0.01252 


0.00894 


-0.00412 


Pz 


0.000422 


0.000301 


-0.000139 


Pz (BP) 


0.000436 


0.000315 


-0.000143 



erties. As in the LDA calculations, we find larger dy 
components and smaller dz components, with a pattern 
of signs determined by symmetry requirements. The in- 
terpretation of the vibron infrared activities in terms of 
this picture has already been discussed in Sec. 0. 

Unfortunately, the level of agreement between the TB 
and LDA results for the molecular dipole moments is 
somewhat disappointing. We find a TB dy value that 
has the right sign, and the correct order of magnitude, 
relative to the LDA value, but the actual values differ 
by a factor of about two. The TB dz value even has the 
wrong sign, but this is related to the fact that the LDA dz 
value happens to come out very small (~10 times smaller 
than for dy-, see also Fig. |6|). Thus, it is not surprising 
that the relative TB error in dz is large, even though 
the absolute TB error is actually smaller for dz than for 
dy. Because the simple TB theory does not include any 
charge self-consistency, it was hoped that the extension 
to the TB+Q theory might improve the results by incor- 
porating a leading (quadrupolar) Coulomb contribution. 
Table Vl shows that the changes from TB to TB-I-Q are 
in the right direction, and there is some improvement in 
the dz (and therefore Pz) values, but the dy discrepancy 
is hardly affected. 

Thus, while the TB theory gives a picture that is quali- 
tatively correct, it is clear that there is room for improve- 
ment. Possible avenues for future investigation may be to 
consider nonorthonormal TB basis sets, to parameterize 
and include off-diagonal matrix elements of the position 
operators between TB basis functions, to include s* or 
d orbitals in the TB basis, or to include other Coulomb 
contributions (e.g., a self-consistent inclusion of the field 
arising from the induced dipoles). 

Finally, we calculated the shell-by-shell spatial decom- 
position of the molecular dipole moments in the TB 
scheme. The results appear in Fig. 04. (By symmetry, all 
dx contributions are identically zero?) Comparing Fig. |lj 
with Fig. |l^, one sees that although the contribution up 
to the first neighbor shell is 91% for the density, it is 



only 62% and 60% for dy and dz , respectively. Up to the 
second shell, the contribution is 97% for the density and 
95% for dy, but only 67% for dz- (Of course, dz is the 
only component that survives in the summation giving 
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removing the orthogonality constraint, which gives rise 
to the long-range tails; the resulting orbitals correspond 
to the "bond orbitals" in Harrison's picture. However, 
the straightforward connection to the polarization is then 
lost, since Pmac is no longer simply the sum of the dipole 
moments of those orbitals, and additional cross-terms 
connecting orbitals in different cells appear. Therefore, 
as far as dielectric properties are concerned, the "maxloc" 
WFs seem to provided the most localized description of 
the electronic structure. In particular, the spatial ex- 
tent of their dipole moments (Figs. ^ and provides a 
natural length scale to compare with the intcrmolccular 
distance in order to assess the validity of the Clausius- 
Mossotti approximation of nonoverlapping dipoles. 



B. Identifying "molecules" in the dense solid 



FIG. 14. Shell-by-shell contribution to dipole. The hori- 
zontal line ( O) is the dipole along the a;-direction. Q s-nd A 
indicate the contributions to dipole along y and z directions, 
respectively. The dipole is in a.u. 



the bulk polarization.) Thus, the dipoles are found to be 
very delocalized, spanning over quite a few neighboring 
molecules, in agreement with the LDA results. 

To summarize this section, we have demonstrated that 
an empirical TB framework allows for a very useful qual- 
itative (and often semiquantitative) analysis of the WFs 
in systems such as the H2 phases under study here. It 
is typical of the empirical TB approach that one can- 
not insist on quantitative accuracy at the level of first- 
principles schemes. However, the TB approximation has 
proven enormously useful over the years because of its 
simplicity, transparency, and ease of application. These 
features often allow for insightful modeling of simple sys- 
tems, or for efficient calculations of large and complex 
systems where ab-initio schemes would not be practical. 
For example, one could easily use the present scheme for 
a computationally efficient analysis of the local dielectric 
structure of more complex H2 crystal structures 
supercell realizations of disordered H2 systems 
expect that the coupling of Wannier and TB approaches 
will prove to be a useful strategy in a wide variety of 
other materials systems as well. 



VIII. DISCUSSION 

A. Wannier functions and the Clausius-Mossotti 
approximation 




There is a vast literature dealing with useful ways of 
identifying individual "atoms" inside a molecule, or indi- 
vidual "molecules" in a dense medium (see Refs. fsj-js^] 
and references cited therein). In this paper we have ad- 
vocated using "maxloc" WFs as a useful computational 
definition of "H2 molecules" in the insulating molecu- 
lar solid (similarly, in Refs. WFs were used for 
defining "water molecules" in liquid water). In spite of 
having some counterintuitive features (becoming larger 
under pressure) and some conceptual limitations (being 
defined only in the independent-electron framework; not 
being totally unique, since they depend on the measure 
of localization), they have the following important con- 
ceptual advantages, (i) When the Clausius-Mossotti ap- 
proximation breaks down - which is when the ambiguity 
in identifying individual "molecules" appears - the sum 
of the dipoles of the Wannier molecules still gives the 
bulk polarization exactly (Eq. ^). By contrast, any def- 
inition of "molecules" based upon a direct partition of 
the bulk electronic charge density necessarily yields an 
incorrect result, since away from the Clausius-Mossotti 
limit the information about the bulk P,nac is not in p{r) 
p3| . For example, in Ref. it was found that differ- 
ent schemes for partitioning the charge density yield very 
different molecular dipoles in liquid water. We would ex- 
pect the situation with such approaches to be even more 
severe in the case of solid hydrogen, since the dipoles are 
smaller and originate mostly in the Wannier tails, (ii) 
At high densities the WFs interpenetrate one another, 
so that the charge density at a given point is a sum of 
contributions from different molecules. Thus, effects re- 
lated to molecular overlap are naturally discussed in the 
Wannier representation. 



As discussed by Harrison [Q, WFs provide a rigor- 
ous formulation of the "extended bond orbitals" which 
appear in tight-binding models. It should be noted that 
orbitals more localized than WFs can be constructed by 
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C. Intramolecular versus intermolecular charge 
transfer 

There has been some debate about whether the strong 
vibron IR activity in compressed sohd hydrogen is due 
mainly to intramolecular or intermolecular charge trans- 
fer |l^j2^,^. As a result of the ambiguity in defining 
"molecules" in the solid, the question is to some extent 
ill-posed. Of course, there is no charge transfer between 
"Wannier molecules", since their charge is fixed. How- 
ever, a heuristic argument can be attempted: the fraction 
of the vibron-induced molecular effective charge originat- 
ing in the central part of the WF can be viewed as the "in- 
tramolecular" contribution (polarization of the molecular 
bond) , whereas the contributions from the orthogonality 
tails in the overlap regions are of "intermolecular" ori- 
gin. As the pressure increases, the "intramolecular" part 
of the WF becomes smaller and contains less charge (and 
presumably becomes less polarizable) , whereas the oppo- 
site happ ens to the outer corona. The results of Sees. 
and VII suggest that as far as polarization-and hence IR 
activity-are concerned, at megabar pressures the "inter- 
molecular" contribution associated with the outer corona 
is dominant. 



D. Summary 

Using a method for computing well-localized Wannier 
functions we have presented a "chemical- like" local- 
ized picture of the electronic structure of solid molec- 
ular hydrogen, and used it to investigate the dielectric 
properties of the compressed system. This approach is 
particularly well-suited for studying the effects of molec- 
ular overlap, which become increasingly more important 
at high pressures. We found the somewhat surprising 
result that already at moderate pressures the orthogo- 
nality tails of the WFs in the overlap regions give rise 
to most of the induced dipole moments on the "Wan- 
nier molecules" ; this clearly indicates a breakdown of 
the Clausius-Mossotti approximation. Under those cir- 
cumstances the electric polarization cannot be extracted 
from the electronic charge density in the unit cell, and 
the Berry-phase/ Wannier- function theory must be 
used instead. The present approach clarifies the origin of 
the strong vibron IR activity in phase III and identifies 
the dominant mechanism: even though the permanent 
dipoles of the molecules in our prototype structures are 
too small to account for the vibron oscillator strength, 
the vibron-induced dipole fluctuations are of the right 
order of magnitude in Cmc2i, and actually too large 
in C2/m. In other words, in the strongly compressed 
solid the dynamical contribution to the vibron effective 
charge dominates the static one. This conclusion seems 
to be supported by the fact that, even though several 
libron modes have been identified in phase III p5| , ^ , 
no strong libron IR activity has been reported. If the 



molecules had any significant spontaneous polarization, 
it should manifest itself in the librational IR absorption. 
Thus, we see that well-localized Wannier functions pro- 
vide a useful definition of "II2 molecules" in the dense 
solid, which can be used to gain important insight into 
the microscopic mechanisms of its dielectric response. 
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